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1. Introduction and General Status 

The main goals of the research under this grant consist of the development of 
mathematical tools and measurement of transport properties necessary for high fidelity modelling 
of crystal growth from the melt and solution, in particular for the Bridgman-Stockbarger growth 
of mercury cadmium telluride (MCT) and the solution growth of triglycine sulphate (TGS). The 
tasks, desribed in detail in the original proposal, include: 

development of a spectral code for moving boundary problems, 

kinematic viscosity measurements on liquid MCT at temperatures close to the (composition 

dependent) melting point, and 

diffusivity measurements on concentrated and supersaturated TGS solutions. 

The work performed during the fifth six-monthly period of this grant has closely followed the 
schedule of tasks. 

The code development has progressed well, though some difficulties have been 
encountered during this period, which are related to the rate of convergence of our spectral 
method. Efforts to remedy this have begun. 

The oscillating-vessel technique that we developed for the measurement of the kinematic 
viscosity of melts at high temperatures and pressures, has been tested with water and gallium. 
Excellent agreement with literature values was obtained. Hence, the technique will be applied to 
MCT in the final reporting period. 

For diffusivity measurements in (aqueous) solutions we have developed and tested a 
technique that overcomes some of the limitations of earlier techniques. During this report period 
we have performed extensive tests with concentrated NaCl solutions. Good agreement with 
published diffusivity data was obtained. Measurements with TGS will be performed in the final 
reporting period. 

In the following we will give detailed descriptions of the work performed for these tasks, 
together with a summary of the resulting publications and presentations. 
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2. MCT Viscosity Measurements. 

A detailed evaluation of the capillary type viscometer for measurements of HgCdTe melts 
was performed during the previous reporting period. Two principal difficulties have been 
identified: 

1. Larger than expected quantities of material (>500g) are needed in order to build up a 
hydrostatic pressure head sufficient to suppress unwanted surface tension effects, such as 
the formation of bubbles in the capillary; 

2. The assessment of safe cell dimensions for high pressure operation is difficult due to the 
complexity of the glassware used in this technique. 

2.1. New Experimental Technique : Oscillating-cup Viscometry 

To circumvent the above shortcomings a complementary technique based on the 
oscillation principle was to be developed. The important advantage of this technique is the 
possibility to use a simple sealed silica cell with a small quantity of material. Two such sealed 
ampules with Hgo. 8 Cdo. 2 T were supplied by Texas Instruments. Hence, we do not have to 
become involved in the preparation of HgCdTe preparation, which would require appropriate 
facilities. In addition, the safety margin with HgCdTe melts in thick-walled silica ampules is 
considerably higher, than in the corresponding capillary method. However, as disadvantages of 
the oscillating-cell or -cup technique one should mention the complexity of instrumentation and 
data analysis. 

In this reporting period, an oscillating-cup or -cell viscometer suitable for viscosity 
measurements of HgCdTe melts has been designed, fabricated, assembled and tested. 

The principle of the oscillating-cup method is based on the damping effect of a contained 
liquid on the oscillatory motion of the container. The energy of the torsional pendulum is 
successively dissipating due to the viscous acceleration and deceleration of the liquid. Thus, the 
time dependence of the angular deflection can be written as 

oc(t) = A exp(- cos(cot + <p) , 

where the angular frequency co = 2jt/T, with T the period of the damped oscillation. The 
viscosity of the liquid can be calculated from the measured logarithmic decrement 8 and the 
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frequency of the oscillation. In addition, the dimensions and density of the liquid, the moment of 
inertia and the logarithmic decrement of the empty suspended system are required for the 
evaluation 
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The most reliable and practical equation for viscosity determination with cylindrical 
oscillating-cup viscometers is the formula derived by Roscoe (1958) 


where 



P 


7tp\ 1/2 

ff, 


R 


30= 1 .1a . -1a 2 
^28 

a^ — 1 + A 1 A 
2 8 


A = 


_s_ 
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The parameters in these relations are as follows: I is the moment of inertia of the suspended 
empty oscillating system, p is the density of the liquid, T is the period of oscillations, R is the 
inner radius of the cylindrical cell, H is the height of the liquid in the cell, ti is the viscosity of 
the liquid and 5 is the logarithmic decrement of the empty torsional pendulum. The above set of 
equations has to be solved numerically for the viscosity. Means of determining the moment of 
inertia I will be described below. 


A schematic view of the developed viscometer is presented in Fig.l. Basically, it consists 
of a torsional oscillation system enclosed in vacuum. A hardened steel (piano) wire of 0.01" 
diameter is suspended from a strain gauge bridge, which, in turn, is suspended from a shaft that 
is inserted through a vacuum feedthrough from the outside. A hand knob at the end of this shaft 
permits manual induction of an oscillation. At the bottom end, the wire is attached to a flywheel 
assembley. This assembley consists of an inertia disk with exchangeable rings. A silica 
suspension connects the flywheel to a sealed silica cell for the liquid. A vacuum system 
(turbomolecular and rotary pumps) can provide a vacuum of the order of 10"5 Torr in order to 
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eliminate damping of the oscillations due to friction with air. The furnace used can be controlled 
with 0.1 °C resolution from "ambient" to 1300 °C. Three independently controlled heater zones 
allow for the adjustment of temperature uniformity of ± 2°C over the length of the silica cell. 

Optical techniques are commonly used for measurements of oscillations. In a standard 
setup, a laser beam, deflected from a mirror attached to the oscillating assembley is projected 
onto a semicircular scale. The amplitude of the oscillation can be determined visually or, 
alternatively, recorded with a photoconductive linear displacement detector. Instead, we are 
recording the torsional amplitude with a simple strain gauge arrangement, and collect and 
process the data in a PC. 

Figure 2 depicts the torsional motion sensing arrangement. Four strain gauges (Omega 
HBM 6/350 LY 11) are glued to the four top stainless steel foil springs of the cross-patterned 
metal strip basket. The torsional deformation of the wire over the length of the basket is coupled 
to the foil springs through the four horizontal spokes of the basket. During the oscillations, these 
foil springs are slightly bent, resulting in a signal from the strain gauges in a bridge arrangement 
with a strain gauge amplifier (Omega model DMD). This particular arrangement is rather 
selective with respect to torsional oscillations and suppresses signals from horizontal oscillations. 

In order to determine a moment of inertia I of the suspended system, the torque constant 
G of the suspended system and its period of oscillation T needs to be measured. Then I can be 
calculated from the relation 

i=g (£F= g “' 2 - 

In order determine G and its possible dependence on the suspended weight, we have performed 
oscillation experiments with 2 different inertial rings added to the flywheel: 

Dimensions and mass of the first ring: R 1=34.42 mm, R2=26. 60 mm, m= 152.4 g 

Dimensions and mass of the second ring: Ri=34.30 mm, R2=26.65 mm, m=261.0 g. 

From the relation for the moment of inertia of a ring about its symmetry axis 

I=m/2(Ri2+R 2 2) , 


we calculated 


1 1 = 1442.2 gcm2 and 12=2462.16 gcm2- 


The measured angular frequencies were: 
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Measured frequencies: 

no rings: wo=1.77 s' 1 > 

with the first ring attached: wi=0.866 s* 1 , 

with the second ring attached: w2=0.699 s _1 • 

From Io = Ii (Oi 2 /((Oo 2 - (Oi 2 ) and Io = h c»> 2 2 /(coo 2 - CO 2 2 ), we obtained Io = 454 gem 2 . 
Furthermore, the three independently calculated values for the torpue constant were: 

G = 454*1.77 2 =1422.3gcm 2 /sec2, G=(454+1442.2)*0.866 2 =1422.06gcm 2 /sec2 and 
G=(454+2462. 16)*0.699 2 =1424.8gcm 2 /sec2. 

Since these values agree well with each other, we can assume that load effects can be ignored in 
the viscosity evaluation in our apparatus. 

The sampling rate of the data aquisition system (DAS) was chosen to be 10/sec, which is 
sufficient for accurate recording of oscillations with periods in the range of 5-10 sec. Total data 
acquisition times of approximately 1 hour were used. The smallest oscillations which our 
detection system can still resolve is about 1 angular degree. Maximum angular amplitudes used 
were of the order of 90 degrees. We have carefully evaluated our data for possible nonlinear 
effects introduced by the detection system or the suspension wire. However, all data consistently 
showed strictly exponential decay of the amplitudes over more than two orders of magnitude; for 
details see below. In addition, the signals were quite sinusoidal with no detectable distorsions 
due to additional parasitic harmonics. The power spectra of the recorded signals showed only a 
single peak and a noise background of a constant magnitude. One third of the total signal trace of 
a typical experiment is presented in Fig.3. 

Initially, the data were analyzed with a routine that calculated the local amplitudes of the 
oscillation as the difference between a local minimum and an adjacent local maximum. 
However, the resulting logarithmic plots of the decaying amplitudes showed often a deviation 
from a straight line in particular at lower amplitudes; see curve 1 in Fig. 4 obtained with water 
(see below). This could be indicative of nonlinear torsional movement, which, in fact, has been 
reported in the literature. Another possible cause is a low signal to noise ratio towards the end of 
the damped oscillation, were linear, pendulous oscillation contributions become relatively 
important. In order to filter the noise, a more complex analysis was performed by fitting a portion 
of the data (usually a couple of oscillations) to a sinusoidal function (nonlinear fitting with 4 
parameters, ao, ai, co, <j>: ao +ai sin(cot+4>)). Amplitudes ai found this way are then plotted on a 
logarithmic scale. They exhibit no marked or systematic deviation from a straight line, as 
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demonstrated by curve 2 in Fig. 4. Consequently, we concluded that no apparent nonlinear 
effects are present under our experimental conditions. 

22. Low Temperature Measurements : Viscosity of Water 

We have performed test measurement with water in a cell in air. The data were processed 
with Roscoe’s formula. The following set of parameters were used : p = 1 g/cm^, R = 1.1 cm, H 
= 16.5 cm, I = 1422*(l/0.9)2 = 1755.6 gcm2, T = 6.981 s, 8=0.0145. These resulted in 
r|=0.00943 poise at 22.5 C, which is in excellent agreement with published data . 

2.3. High Temperature Measurements: Viscosity of Molten GaUium. 

A test experiment with molten gallium was performed to evaluate the heating 
arrangement. First, 100.05 g of gallium in form of pellets was loaded into the glass cell. Then, 
the logarithmic decrement was measured in air as well as in vacuum. The following values were 
obtained: in air: 8=0.00331; in vacuum 8=0.00222. After taking these measurements, the 

furnace was powered. In Fig. 5 one can clearly discern the moment when gallium began to melt. 
Measurements were made at various temperatures between 27 and 300 °C. Figures 6-8 show the 
plots of amlitudes versus time obtained at 200, 250 and 300 °C, respectively. The measurement 
at 200 °C was repeated once after a complete temperature cycle and gave the same result within 
1% . For the viscosity calculations of gallium, the following data were used: (0 = 0.855 sec'l> p 
= 6.1 g/cm^, R = 1.1 cm, H = 4. 2cm, G = 1422 gcm^/sec^. The 8’s obtained at the various 
temperatures, after subtraction of the intrinsic logarithmic decrement of the system (8 va c) are 
listed in Table 1. 


Table 1: Logarithmic decrements versus temperature obtained for gallium 
Temperature PCI 


27 

0.01406 


0.01227 


0.01147 

200 

0.01125 


0.01082 


HXirem 


The resulting viscosities are plotted in Fig. 9, together with literature data which, assuming 
Arrhenius behavior, follow the relation Tj = 0.3567*10217. 4/T C p While there is good 

agreement below 150 °C, there is a deviation of about 10% at the higher temperatures. 



In conclusion, the equipment for viscosity measurements of molten HgCdTe was tested 
using water and liquid gallium. The overall performance is quite satisfactory. We expect to 
obtain viscosity data on HgCdTe system with an accuracy of at least 95%; the temperature slope 
of the viscosity, which is particularly important for theoretical evaluations can be expected to be 
resolved with even higher accuracy. 


3. TGS Diffusivity Measurements 

In earlier reports our approach to measuring the diffusion coefficients of supersaturated 
solutions was outlined. A rectangular optical cell is initially partly filled with a solution of a 
certain concentration C. A solution of higher concentration C + AC is then injected at the bottom 
of the cell with a syringe. The resulting system is convectively stable with the heavier solution 
below the lighter one and mixing occurs by diffusion only. A Zygo Mark II Mach-Zehnder 
interferometer interfaced with a personal computer is used to follow the evolution of the 
concentration profile in the cell. At regular intervals the interferometric intensity profiles 
produced are stored. The advanced fringe analysis software ZAPPC is then used to convert these 
into refractive index profiles which are proportional to the concentration. A numerical 
integration of these profiles, described in detail in the last progress report, yields the diffusivity. 

In order to assume that the diffusivity D is a constant in the range of AC and keep the 
evaluation method numerically tractable, this concentration difference has to be small. 
Additionally, the supersaturated concentration region extends from the saturation point to the 
spinodal limit and this range is small for most solutions, usually well below 1M. For 
concentration measurements in this region AC has to be smaller than this range. These 
considerations suggest that it is preferable to maintain AC at 0.1M. 

The effect that runs counter to this is the convection which is caused by the injection of 
the solution of concentration C+AC into the cell. When AC is large the density stratification in 
the cell is more stable and convection caused by injection subsides faster, producing less scatter 
in the data. Hence, in order to obtain good reproducibility the two competing effects had to be 
balanced. With a AC of 0.1 M convection persists for longer times as the heavier solution now 
settles more slowly, with the resulting scatter in the data as high as 5%. To dampen convection 
without resorting to excessive values of AC, progressively thinner cells were employed. Our 
criterion was that the data should have scatter of less than 2% for a AC of 0.2M. In the last 
report results were given for cells of 1mm light path and they were still somewhat unsatisfactory. 
In the last few months new cells of smaller path lengths were used in an attempt to improve these 
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results further. Cells of 0.2mm path length completely damped out convection but now the path 
length was so small that it was difficult to observe and measure the concentration differences 
interferometrically, that is the effective optical path lengths differences between various points in 
the cell were very small. Cells of 0.5mm path length gave a good balance between measurability 
and dampening of convection. 

Unfortunately cells of less than 1mm path length are not manufactured in standard 
rectangular shape. Additionally, due to the difficulty of cleaning cells that thin, they are made 
with detachable covers and in our experiment this led to problems of seepage between the glass 
plates. Cells of 1mm path length inserted with microscope slides to reduce the path length were 
also tried and here again it was difficult to avoid seepage of the solution between the slides. All 
this meant that the process of obtaining good data had become very difficult and time 
consuming. At this point the decision was made to abandon the optical cells provided by cell 
manufacturers and design our own system. 

3.1 New Interferometry Cell 

An exploded presentation of the new demountable and readily cleanable cell is given in 
Fig. 10. The cell consists of two rectangular windows that sandwich a Viton gasket. Four bolts 
excert sufficient pressure onto the windows through two precision-machined brass flanges with 
highly planar faces to prevent any leakage. The windows consist of standard microscope slides 
that have been tested interferometrically. The solutions are injected through the open top and 
viewed through the windows in the flanges. The light path length is equal to the thickness of the 
rubber sheet and can be adjusted by using sheets of different thickness. The results reported 
below were obtained using rubber of 0.7mm thickness. 

With this arrangement all the difficulties associated with this method, as reported earlier, 
have been overcome. Diffusivities of supersaturated solutions can now be measured with good 
reproducibility. An additional benefit is that only flat glass slides are used here, not glass cells 
which a have greater probability of having stress points. This lack of stress points in the glass 
produces fringe patterns of greater uniformity. All this is shown by the reproducibility of the 
data which is better than ±1.5%. 

3.2 High Concentration NaCl Solutions 

Before application of this technique to TGS, we have performed careful measurements 
with saturated and supersaturated solutions of NaCl, for which data are available ion the 
literature. Our results are summarized in Table 1. Note that the saturation point at 25 °C is 5.3 



molar (M). The standard deviations shown are the result of six independent measurements for 
each data point. 
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Table 1. Measured Diffusivities of Concentrated NaCl solutions at 25 °C 


Method Used 


Diffusivities (cm 2 /s) 



5M 

5.2M 

5.4M 

Previous work [1, 2 & 3] 

1. 584xl0- 4 5 

1. 580xl0' 5 

1.49x10-5 

Light path=lmm, AC=0.15M 

(1.56±0.04)xl0- 5 



Light path=lmm, AC=0.3M 

(1.57±0.02)xl0’ 5 



Light path=0.7mm, AC=0.2M 

(1.58±0.02)xl0- 5 

(1.57±0.02)xl0-5 

(1.48±0.02)xl0-5 


At this stage we can state that an accurate new interferometric method of measuring 
diffusion coefficients of saturated and supersaturated solutions has been developed. The only 
other attempt to measure diffusivities of supersaturated solutions was made by Myerson and co- 
workers [3]. Our method is not only simpler but provides better reproducibility than the ±3% 
achieved in [3]. Work on the diffusivities of supersaturated TGS solutions utilizing this 
technique will begin shortly. 

References 

[1] V. Vitagliano and P.A. Lyons, Diffusion Coefficients for Aqueous Solutions of Sodium 
Chloride and Barium Chloride, J. Amer. Chem. Soc., 78 (1956) 1549. 

[2] J.A. Rard and D.G. Miller, The Mutual Diffusion Coefficients ofNaCl-H20 and CaCl 2 - 
H 2 O at 25°C from Rayleigh Interferometry, J. Solution Chem, 8(1979) 701. 

[3] Y.C. Chang and A.S. Myerson, The Diffusivity of Potassium Chloride and Sodium 
Chloride in Concentrated, Saturated and Supersaturated Solutions, AIChE J., 31 (1985) 
890. 

4. MCT Code Development 
4.1 Introduction 

During the last semi-annual period work was completed on the solution of the moving 
boundary problem for axisymmetric solidification in the absence of convection. In addition, a 
code which solves for axisymmetric convective flow in a circular cylinder using a Chebyshev 
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pseudo-spectral method together with an influence matrix method was also completed. The 
solidification code development is described in section 4.2 and the axisymmetric convection code 
is described in section 4.3. 

4.2 Solidification: Code Development 
Governing equations and boundary conditions 

The governing equations and boundary conditions are described in terms of a cylindrical 
coordinate system with the orgin located at the centerline and the top of the ampoule. Variables 
are put in dimensionless form by scaling lengths with the ampoule radius R and defining the 
dimensionless temperature as T=(T - T c )/(Th - T c ) where Th - T c is the (fixed) temperature 
difference between the top and bottom of the ampoule 


AT m =Pej|^, 

(1) 

AT c =Pe c |^, 

(2) 

AT a =Pej^, 

(3) 


where subscripts m, c, and a refer to melt, crystal and ampoule, respectively. Here, the Peclet 
number, Pe =V a R/a, represents the dimensionless transition rate of the ampoule, and V a , R and 
a are, respectively, the growth velocity, ampoule radius and thermal diffusivity. 

Boundary conditions 


T(r, 0) = 1 , 

(4) 

T(r, A) = 0 , 

(5) 

3T(0,z) 
dr ~ U ’ 

(6) 

aT(Rw,z) =Bi(TTdi 

(7) 
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where A = L/R is the aspect ratio of the ampoule, with L being the ampoule length, R w = (R + 
d w )/R with d w = ampoule wall thickness, Bi(z) is a dimensionless heat transfer coefficient 
defined to include radiative and convective heat transfer between the ampoule and the furnace, 
and Tf (z) is the temperature distribution of the furnace wall. 

At the melt/solid interface 


9T 

9n 


Krr 

K c 


-)T 

3 - = Ste Pe(e z - n) , 

dn m 


( 8 ) 


v T = T* , (9) 

where n is the unit vector normal to the interface, e z is the unit vector in cylindrical coordinate 
and Ste = AH/Cp AT is the Stefan number and T* is the dimensionless melting temperature of the 
crystal. Perfect thermal contact at the melt/ampoule and crystal/ampoule interface is assumed. 
Thus, the heat fluxes at these boundaries are equal i.e. 


im .M3 and(a=§©, 
l dr /m K m \ dr ; a \ or ) c Kc\ dr / a 


( 10 ) 


Solution method 
Nonorthogonal transformation 

The nonorthogonal coordinate transformations utilized in this study (Fig. 11) are the 
following: 


Upper domain, % = r, T[=-f - , 

h(r) 

Lower domain, \ = r, rj=2 — z ) 

h(r) - A 


where h(r) is the melt/crystal interface position. In the ampoule, h(r)=h(R)=constant. The 
transformed governing equations now take the form 


9 2 T .^T 

+ A 

dZ? ^i 2 








(ID 


where in the upper domain, 
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A=j- + (a^f 


h 2 l h 


B = - 


2r|3h 

h 3r 


C = — , 


lh dr) h-\ 2 tu dr h 


l 3r z £h 

Here, Pe should be taken as Pe m and Pe a , respectively, for the melt and ampoule. 


In the lower domain, 


A = — 1 — +1 

'1-2 3h\ 2 

M* 1 

^h-A ^1 ’ 

D _ 2(n-2)3h 

c = 1 

h-A 3r 

$ 

f i ah) 2 i a 2 h 

i ah, 

ih-A dr l h-A3r 2 

£(h-A) 


where Pe should be taken as Pec and Pe a respectively in the crystal and ampoule. The boundary 
conditions are also subject to the coordinate transformations. 


Pseudospectral discretization 

The total solution domain is divided into four subdomains (Fig.l) and for each subdomain 
the equation (11) is discretized using a Chebyshev-collocation method. The Chebyshev 
polynomial expansion for the dependent variables is 


Mjsr 

T =I cuWtCy)), 


i=l 

j=l 


( 12 ) 


where tj(x) = cos(iarccos(x)), tj(y) = cos(j arccos y), and x and y are defined by the 
transformations given in Fig. 1 1 for each of the subdomains. The Gauss-Lobatto collocation 
points are 


x=cos ' / 
1 M-l 


and y.=cos-^— ), 


N-l 

The Chebyshev collocation derivative of a function T(x,y) can be represented as 


( 13 ) 
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(§R -2 Dj’rCXp.yj), 

' ox 'y p=i 

(f) 

' y, V q=l 


y 1 J q=l 


[dY 


\Py 2 1 


= I ^(x,, yj , 


q=l 


->2_ \ m,n 

3 T ' - 2 ^EfT(x p ,y q ) • 


\dxdy).. p=! 


'y 


q=l 


( 14 ) 

(15) 
(15) 

(17) 

(18) 


where the expressions for D x , D y , D xx , D yy and D xy are given by Ouazzani[l].For large m, n, 
the Chebyshev collocation differentiation can be efficiently implemented by using FFT’s (Fast 
Fourier Transforms), which eliminates the use of expensive matrix multiplication involved with 
(14)-(18). 


PGCR method for finite difference preconditioning 

In the present investigation, the preconditioning solution method for spectral equations 
involves finite difference preconditioning imbedded within each iteration. This allows the 
problem to retain spectral accuracy while using a finite-difference method to solve a related 
auxiliary problem. For this purpose, we need a finite difference discretization for the 
preconditioned version of equation (11). The following second order accurate finite difference 
formulas are employed for the entire non-uniform grid: 

W - dX; fj-M-fj , dx i+ i fj-fj,! 

[dx/i dxj+dxi+i dxi+i dxi+dxi+i dx; 


dV 

’^x 2 , 


= 2 


i+1 


i-l 


f. 


[ (dx.+dx i+1 )dx i+1 (dx+dx +1 )dx dx.dx 


(19) 


1 + 1 / 
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The resulting finite-difference equation takes the form: 

Acij*t*ij + Aejj^i+ij +Aw;j(j)j.i i j -t-Anij^jj.).] +As;j +Aenjj(j)j + |j + i +Aesjj<{)j + i i j_i 

+Awn..<ti , . . +Aws..<jx , . , =R. (20) 

where Rjj is the residual of equation (11) at node (i,j) and (J)!! is the difference in the variable T 
evaluated at the current and previous steps. This leads to 

A <j> = R. (21) 

The finite-difference matrix A is a nine-diagonal coefficient matrix. It is approximated by 
incomplete factorization i.e. A=L U , where L is a lower triangular matrix and U is an upper 
Triangular matrix, and the factors are easy to calculate. The approximate factors L and U may 
have the sparsity structure of the original matrix. The solution for equation (21) is 
straightforward since it consists of forward substitution 

U*<t> = (L*)-lR, (22) 

followed by a backward substitution, - 

4> = (UV (L*)-*R. (23) 

The subsequent approximations are obtained via 

T m+1 =T n +X(J) , (24) 

where X is a relaxation parameter. Unfortunately, such an iteration scheme generally converges 
slowly, so an acceleration technique must be added. Alternative acceleration techniques which 
have been developed are conjugate direction methods. For sparse symmetric positive-definite 
linear systems the preconditioned conjugate residual method is one of the best. However, 
equation (11) is far from symmetric, therefore, in the present investigation the preconditioned 
generalized conjugate residual (PGCR) method is the chosen acceleration technique since it was 
developed specifically for nonsymmetric systems. 

The PGCR procedure is as follows: 

Initialize 

To 

R°=F - LT° 

A<}> 0 = R° 

P° = <() 0 



Iterate 
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_ (R m , LP m ) 
(LP m ,LP m ) 


T m+1 = T m + oc m+1 P m 


R m+ i = R m - o^ +1 LP™ • (25) 

A(J) m+1 = R m+1 

_ (L(j) m+1 ,LPi) 

J (LPi,LPj) 


pm+l _ ( j ) m+1 + £ pj“ +1 pj 
j=0 


The scalar a is chosen to minimize the Euclidean norm of the residual. The P's are 
orthogonality coefficients and result from the orthogonality property for the conjugate residual 
version. The practical calculation only requires the following condition to terminate the iteration 


max 


|'pn+l _ 

'pn+1 




(26) 


where we took e = 10~7 • 

The PGCR method is an attractive iterative method for solving non-symmetric system 
since the iteration cannot break down, it is stable and the minimization increases the convergence 
rate considerably. We compared the PGCR method to a Conjugate Gradient Squared (CGS) 
method and a Strongly Implicit Procedure (SIP). For all these methods, we also used incomplete 
factorization to solve the preconditioned problem. The results show that PGCR method is more 
efficient in terms of CPU time. 

Outer iteration 

Having computed the new values of temperature, an outer-loop iteration is required to 
update the melt/solid interface h(r) at which the computed temperature is equal to Tmelt. One 
scheme to update the interface shape is Newton's method 


h, L+1 




'i 


(27) 
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where Fi = Ty* - T me i t , here Tjj* is the calculated temperature at the Lth interface location. 
However for Newton's method, it is difficult to guarantee the outer iteration convergence for all 
cases. Therefore as an alternative to Newton's method, an interpolation scheme is used to update 
the interface. In the present study, the new interface is located by the following linear 
interpolation: 

if T m < Tjj and T m > T ; >j+1 

-Zj). (28) 

1 1 ij ‘ 1 i.j+1 ) 


The outer-loop iteration is continued until the following convergence criterion is satisfied: 


where we took e = 10" 4 . 


ER = ma» 


h L+1 - h L 


h L+1 




(29) 


Results and discussion 

For the purposes of comparison of our method with that of Adomato and Brown [2] (for 
cases in which convection did not influence the interface shape) we have carried out numerical 
simulations for the growth (in the absence of convection) of a gallium-doped germanium (Ga:Ge) 
alloy in a vertical Bridgman system. (Thus, our computed interface shapes, which are for 
solidification in the absence of convection could be compared with theirs). The calculations for 
Ga:Ge growth are presented separately for the two furnaces and the different ampoule. For 
a heat-pipe furnace, the changes in the heat transfer coefficient between ampoule and three 
zones of furnace are modelled by the function 

Bi(z) =^p(c 3 {l + tanh[ci(zc - c 2 - z)]}+l+ tanh[-ci(z c + c 2 -z)]) , 

where zc is the location of the mid-plane of gradient zone. The Bio is the value of Biot number in 
the cold zone, C 3 sets the ratio between the Biot numbers in the hot and cold zone, ci is the slope 
of the transition in Bi (z) between the isothermal zones and the adiabatic region, and c 2 is the 
half of the dimensionless length of the adiabatic zone. The furnace wall temperature is modelled 
by the function 


Tf(z) =0.5{1 + tanh[ 12(0.5 - z/A)]}, 
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where the constants in this expression have been fit to the measured wall temperature of the 
furnace. For details of the thermophysical properties and growth conditions the reader is referred 
to [2]. 

For a heat pipe furnace and a boron nitride ampoule, the temperature field is shown in 
Fig. 12a. The difference between the thermal conductivities of the melt, crystal and ampoule 
cause the melt /crystal interface to be convex with respect to the melt and results in the 
temperature decreasing radially adjacent to the phase boundary. The mismatch in thermal 
boundary conditions between the adiabatic and hot zones of the furnace causes a second set of 
radial gradients with the hottest temperature located at the ampoule wall. Fig. 12b depicts the 
temperature field using a quartz ampoule (in a heat pipe furnace) which decreases the 
conductivity of the ampoule to a factor of 14 lower than the melt. As seen by Adomato and 
Brown, the ampoule wall takes up most of the radial temperature difference required to transport 
heat in the system. Thus, the interface is flatter, and the radial temperature gradients near the 
interface are smaller. Larger radial temperature gradients occur near the intersection of the 
adiabatic and hot (cold) zones. 

The constant gradient furnace is modelled by using a constant Biot number over the entire 
length of the ampoule and specifying the furnace temperature profile as 

T f = 1 - z/A . 

The isotherms are almost flat, the large radial gradients that arise at the intersection of adiabatic 
and hot (cold) regions in the heat-pipe furbace is eleminated, and the deflection of the 
melt/crystal interface is large as indicated seen in Fig. 2.c. 

In addition to the work described above, we have started to add convection in the melt 
phase. The biggest problem encountered to date appears to be that of convergence rates which at 
present are excessively slow even with the finite difference preconditioning. 

Future Work 

Work planned for the next six monthly period includes the incorporation of melt 
convection into the model. Based on our experience with convergence rates for the problem of 
convection in a cylinder (see section 4.3) using the influence matrix method (without 
preconditioning) the prospects of obtaining an efficient solution method for the coupled 
convection-solidification problem are presently not good. We will thus continue to investigate 
convergence acceleration techniques for spectral methods. 
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4.3 The Influence Matrix Method Applied to Convection in a Cylinder 
Governing equations 

For an axisymmetric system, in cylindrical coordinates the heat and momentum transport 
equations in stream function - vorticity form may be written as 


dt 

^ + U r ^ + U z ^i-^L=v(^- + l 
dt dr dz r \ dr 2 r 


o cue i dcoe < 

+ “ d? + a?2 


\2~ 


d y i dy d y 
— — -zr~ + — — = rcoe, 
a? 2 r 3r 3^2 


( 1 ) 

( 2 ) 

( 3 ) 


where T represents temperature, oie the azimuthal component of the vorticity, y the stream 
function, t time, a the thermal diffusivity, v the kinematic viscosity, (3 the coefficient of thermal 
expansion, g z the axial component of the gravitational acceleration, u the velocity vector, and 
u Zj u r the axial and radial components of the velocity, respectively. The " indicates a 
dimensional variable. 

The stream function y in eqn. (3) is defined by 


u r 



r 3z ’ 


(4) 


and the azimuthal component of vorticity is defined by 


~ 3u r 3u z 
^ 3z 3r 

For an axisymmetric system the two other components of vorticity are identically zero. 


(5) 


Due to the extra diffusive component in eqn. (2) c D©/? 2 , which arises in cylindrical 
coordinates, an alternative variable S = rcog must be solved for so that the discretization in time is 
consistent (details given below). Therefore 


dS ~ dS ~ 3S 

-TT + u r3=- + U z ^r ■ 
3t dr dz 


^ u r = v 


f d 2 s 1 dS cfs_ ] 

dr 2 ? d? + a ^2 


dT 


(6) 



20 


and 


dr 2 ? 3t 0^2 


are solved for in place of eqn.s (2) and (3). 


( 7 ) 


Equations (1), (6), and (7) are solved inside a cylinder of radius R and height H. Due to 
axisymmetry, the boundary condition 


9T 

dr 


= 0 , 

r = 0 


( 8 ) 


is applied at the centerline. Dirichlet, Neumann, or Cauchy boundary conditions may be applied 
for the temperature at the top, bottom, and side of the cylinder. 


The no-slip condition for the fluid velocity at the top, bottom, and side of the cylinder 
results in the boundary conditions 


d\\r 

dz 


z=— 


-o* 

dz 


z = 


= 0. ^ 

dr 


= 0, 


r = 0 


(9) 


for the stream function. Since the value of the stream function is arbitrary to within a constant, its 
value is set to zero at the top, bottom, and side of the cylinder. This combined with the condition 
of no penetration for the fluid velocity at the centerline determines the value of the stream 
function at the centerline to be zero. 


Because of axisymmetry, S=0 at the centerline. As is the case for vorticity in the stream 
function - vorticity formulation, the boundary conditions for S at the other three boundaries are 
not predeterminable. 


The general approach of scaling the velocities by a characteristic velocity U, the times by 
a characteristic time t, and lengths by a characteristic length L results in the nondimensional 
equations 


(Fo) 4 ^+PrRe |u, 


d© d© 

1T + Uz d7 


d 0 i d© d 0 
+ T + , 

dr 2 r dz 2 


( 10 ) 


(FoPr)’ 1 ^- + Re[u 


dS 

dr 


„ as 2Su r ) 
Z 5z ‘ ~l 


3 2 S ; 8S ( S 2 S r Ra 80 
0 r 2 r dr 0 Z 2 Re Pr dr 


dV 1 , d 2 y ^ g 

dr 2 r ^ dz 2 


( 11 ) 


( 12 ) 



where Fo, Pr, Re, and Ra are the Fourier Number, the Prandtl Number, the Reynolds Number, 
and the Rayleigh Number. The definitions of these four nondimensional quantities are given in a 
table at the end of this section. The three field variables are defined by 
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(13) 

AT U UL 2 

where AT is a characteristic temperature difference andT re f is a reference temperature. 

In nondimensionalizing the boundary conditions it is convenient to choose the 
characteristic length to be half the radius (i.e. L = R/2). With this choice the boundary conditions 
for equation (10) are 

50 

dr 

and at the top (z = +1/L), bottom (z = -1/L), and side of the cylinder (r = 2) Dirichlet, Neumann, 
or Cauchy conditions may be applied. Note that A = R/H is the aspect ratio. For equation (11) 
the boundary condition at the centerline is S=0, and the boundary conditions at the top, bottom, 
and side of the cylinder are not predeterminable. The boundary conditions for equation (12) are 

= 0, (15) 


9y 

= 0, & 

= o. *£ 

dr 

r= 2 3z 

z = +l/A 


r = 0 


= 0 , 


(14) 


y = 0 at r = 0, y = 0 at r = 2, y = 0 at z = 1/A, y = 0 at z = -1/A, (16) 

Numerical Method 

Equations (10), (1 1), and (12) are discretized in time using a semi-implicit scheme for the 
diffusive terms and an explicit scheme for the convective terms. Rearranging the discretized 
equations into the form necessary for their solution by means of a Helmholtz solver yields 


l d 2 @ 


d 2 ©\ 


1 90 

+ jr - 5 - + 

I9r 2 r ar 9z 2 


,n+l 


2 At Fo 


0 ' 


n+1 


0 nl . 40 n , J 30 , 30 f I 30 , 9©f 4 

LT'3T +u ^) -br +u ^J . 


( 17 ) 
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S n-1 _ 4S n 
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2 S u r ' 

r 

1 as 

as 

r a7 +Uz aT' 

r J 

1 - 

h* +U: 

l dz " 

aV +1 

i a\j/ n+1 

- 4 - 

aV +1 

s n+1 

3r 2 

r ar 

i 

az 2 


S n+1 


l&n + j: 

r | R< 


n+1 


Ra d®' 
Re Pr 5r 


(18) 


(19) 


where the integers n-1, n, n+1 indicate the number of the time step. 


Due to the lack of boundary conditions for equation (18), the method of the matrix of 
influence is used for the solution of the momentum part of the problem, that is equations (18) and 
(19). In the method of the matrix of influence first a preprocessing stage must be completed. In 
this stage a set of elementary solutions are calculated. These elementary solutions are defined by 


d% , as, 1 3 2 s, , : Q 

3r 2 r ^ 9z 2 2 At Fo Pr 

with boundary conditions 

SjOli) = 8ij> 'Hi 6 Tmi, 

SjOli) = Oi Tli £ r- Fj 4 i, 

Sj(rii) = 0, r = 0, 
and 

dVj 1 dVi , 8 Vj _ c 

3r 2 " r ^ 3z 2 j> 

with boundary conditions 

¥j = 0, at r = 0, r = 2, z = -L, z = — , 

A A 


(20) 


( 21 ) 


( 22 ) 

(23) 


where 8ij is the Kronecker delta, T li are the collocation points that lie on the boundary of the 
computational domain, T represents all of the boundary except for the centerline, andrMi is the 
special matrix of influence subset of T. The final two steps of the preprocessing stage are the 
generation of the matrix of influence, and then to take its inverse. The matrix of influence is 



Tli e Tmi, 


(24) 


where n is the outward pointing unit normal at each boundary. Due to the occurence of At in 
equation (20) the preprocessing steps must be repeated if the size of the time step is changed. 


'i 
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At each time step seven stages are completed. The first stage is solution of 


9^* 

as 

a r 2 


ids a^s 

r 3r + 


dz 2 2 At Fo Pr 


_ S"- 1 - 4S n 
2 At Fo Pr 


Re 


j as ^ as 2 s u> | as as 2 s mV 1 - 1 
Tar + U *3F • -T 1 ) - hr + • -r 1 ) . 


n+1 


r Ra 

Re Pr 3r 


(25) 


with boundary conditions 


S = 0, at r = 0, r = 2, z = — , z = — . 

A A 

The second stage is solution of 

d 2 ¥ i ay a 2 y $ 
a r 2 " r ^ a z 2 ’ 

with boundary conditions 

y = 0, at r = 0, r = 2, z = — , z = — . 

A A 

/N 

Using y from the second stage, in the third stage h is generated: 


= ’Hi e F M i • 

3n 


In the fourth stage A, is generated: 


h = H^i) = M- 1 ^ h(rij) . 


The Fifth stage is the solution of equation (18) with boundary conditions 


(26) 

(27) 

(28) 

(29) 

(30) 


S n+ 1 (Tii) = A if Tjj e T mi , 

s n+1 ("Hi) = o, tjj 6 r - r mi* (3i) 

S n+1 (Tli) = 0, r = 0, 

In the sixth stage the stream function at the (n+l)* time step is generated by solving equation 
(19) with boundary conditions 


y n+1 = 0 at r = 0, y n+1 = 0 at r = 2, y n+1 = 0 at z = 1/A, y n+1 = 0 at z = -1/A. (32) 

Finally, in the seventh stage the true values of S n+1 on the boundary are calculated by evaluating 
the left-hand side of equation (19). 
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Non-Dimensional Groups 


Fourier Number 

Fo = a x / L 1 2 

Prandd Number 

Pr = v/a 

Reynolds Number 

Re = U L / v 

Rayleigh Number 

Ranp g ATL 3 /va 
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